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Abstract—In this paper, the periodic analysis of candidate plas- 
monic topologies for the implementation of left-handed media at 
optical frequencies is pursued through a triangular-mesh-based 
finite-difference time-domain (FDTD), equipped with Floquet 
boundary conditions. The technique is shown to possess excellent 
convergence and accuracy properties, as opposed to the conven- 
tional rectangular-cell-based FDTD. The latter fails to accurately 
capture the plasmonic resonant modes excited in the lattices under 
consideration. The studies presented in this paper are particularly 
aimed at rigorously investigating the possibility of backward-wave 
propagation in periodic arrays of plasmonic nanoparticles, along 
with their potential analogy to microwave negative-refractive 
index transmission-line metamaterials. 


Index Terms—Finite-difference methods, left-handed metama- 
terials, plasmonics. 


I. INTRODUCTION 


TIS WIDELY recognized that Veselago’s early work on the 

possibility of media concurrently exhibiting negative 
dielectric-permittivity « and magnetic-permeability jz: [1] has 
been the main inspiration for the recent activity in the area 
of metamaterials. Veselago introduced the term left handed for 
such media after the formation of a left-handed triplet by the 
electric, magnetic, and wave vectors in these and suggested that 
the index of refraction in left-handed media would be negative. 
Thus, a number of unconventional phenomena, including neg- 
ative refraction at a positive to negative index interface, would 
be enabled. 

The theoretical results of [1] were experimentally confirmed 
through microwave artificial dielectrics exhibiting negative- 
refractive index (NRI), with the split-ring resonator and strip- 
wire (SRR/SW) structure of [2] being the first realization of 
such. A planar alternative exhibiting a wide NRI bandwidth 
compared to the strongly resonant SRR/SW geometry was 
proposed in [3] and [4]. The latter was based on a 2-D transmis- 
sion line (TL) loaded by chip or distributed implementations 
of lumped elements. Therefore, the extension of the NRI-TL 
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concept to optical frequencies, while appealing due to its 
potential for broad bandwidth, becomes synonymous with the 
conception of optical analogs of lumped elements. 

Recently, Engheta et al. addressed this question by identify- 
ing lattices of plasmonic spheres/ellipsoids as possible counter- 
parts of the NRI-TL structure in the infrared and visible regime 
(5]-[7]. Furthermore, 1-D and 2-D arrays of silver nanoparti- 
cles have been shown to exhibit backward-wave bands, which 
is the signature property of NRI media, within their Brillouin 
zone [8], [9]. In [9], a 2-D silver-rod lattice of a deeply sub- 
wavelength unit cell was studied. A plasmonic backward-wave 
mode found in the Brillouin zone of this “subwavelength 
photonic crystal” (SPC) structure was employed to achieve a 
realization of Pendry’s “perfect lens” [10]. Interestingly, the 
lens structure of [9] presented a very small bandwidth of about 
0.01%, while similarly intriguing questions of analysis, design, 
and optimization arise as candidate NRI topologies are consid- 
ered at optical frequencies. 

In this paper, the dispersion analysis of periodic structures 
of plasmonic rods is performed through the finite-difference 
time-domain (FDTD) technique [11], [12]. FDTD is simple, 
versatile, and provides for the wideband characterization of 
geometries in a single simulation. However, its Cartesian mesh 
version would approximate a cylindrical or spherical inclusion 
by mesh staircasing. In principle, the well-known numerical 
errors due to staircasing [13], [14] can be reduced by mesh 
refinement of the standard FDTD or by modified FDTD 
schemes aimed at improving the modeling of nonrectangular 
boundaries [15]-[17]. The studies included here show that 
these errors become quite pronounced when plasmonic modes 
are involved, resulting in spurious resonances. In a dispersion 
analysis, these resonances would manifest themselves as spu- 
rious, mesh-dependent, and poorly convergent bands. Similar 
effects have been observed in the application of FDTD to scat- 
tering problems with significant surface-wave contribution [18]. 
To overcome this problem, a formulation of FDTD in a tri- 
angular mesh is adopted that is similar to the simple explicit 
schemes proposed in [18] and [19]. The formulation is extended 
to incorporate periodic boundary conditions, thus reducing 
the computational domain to a single-unit cell of a periodic 
structure via the sine—cosine method [20], [21] and applied to 
characterize plasmonic structures of interest. 

In the following, our approach to the FDTD dispersion 
analysis of 2-D arrays of plasmonic nanoparticles is presented, 
along with a discussion of its numerical stability properties and 
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Fig. 1. Triangular FDTD mesh. 


a comparison with the conventional rectangular-mesh FDTD. 
Subsequently, the proposed methodology is employed for the 
calculation of the complete band diagram of the SPC of [9], 
illuminating the rich band structure of the latter and providing a 
solid background for the understanding of its operation, includ- 
ing its narrowband behavior as a lens. Finally, the possibility 
of the SPC lattice operating as an optical analog of a recently 
proposed microwave negative-refractive-index medium [23] is 
investigated. 


Il. TRIANGULAR-MESH FDTD FORMULATION 
A. Update Equations 


A 2-D transverse-electric case is considered, with E,-, E,-, 
and H,-field components. The computational domain is divided 
in triangular elements (as shown in Fig. 1) and generated 
by a quality Delauny triangular-mesh generator [24]. In each 
element, tangential electric-field components are sampled at 
the center of each edge and the perpendicular magnetic-field 
component is sampled at a “centroid” point, whose coordinates 
are computed by averaging the coordinates of the triangle 
vertices. In [18] and [19], 1, was sampled at the circumcenter 
of each element, which may possibly lie outside a triangle con- 
taining an obtuse angle. In that case, the distance between /7, 
sampling points in adjacent triangles may diminish, aggravating 
the stability of the resulting FDTD scheme. This problem is 
overcome by sampling H, at the centroid instead. Note that 
the two approaches tend to become equivalent in the limit of 
a dense mesh, where the Delauny mesh generator is expected 
to produce nearly equilateral triangles with almost collocated 
centroids and circumcenters. 

Upon setting up the mesh, the derivation of the field update 
equations proceeds as follows. For the update of A, the 
integral form of Faraday’s law is discretized. Considering the 
geometry of Fig. 2, the line integral of the electric field is com- 
puted along the triangular contour enclosing the element shown, 
while the surface integral of the magnetic field is approximated 
by the product of the H, component at the centroid of the 
triangle times its area: 


OH, 
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Fig. 2. Geometry demonstrating the update equation for a perpendicular 
magnetic-field component. 


P, 


Fig. 3. Geometry demonstrating the update equation for a tangential electric- 
field component. 


Approximating the remaining time derivative by a second-order 
accurate centered difference, the following update equation is 
derived: 


— At me im a 
Heth? yn MO + (Ee |PaP al + Be, [PsP + Be, |P.Pa) 
(2) 


where At is the time step, and n is a time-step index. Similarly, 
the update equation for the tangential electric-field components 
stems from the application of Ampere’s law, using the line 
segment connecting the centroids of two adjacent triangles, as 
shown in Fig. 3, as an Amperian path: 


oD: 
Ot 
Hence, the update equation for a tangential electric-flux-density 
component D; assumes the form 


=) _ 
|Q1Q2| sin t 


The incorporation of the material dispersion characterizing 
the plasmonic nanoparticles is performed by translating the 
frequency-domain constitutive relation D;,(w) = e(w)E;(w) 
(where ~ denotes a phasor quantity) to the time-domain. To 
that end, an auxiliary differential-equation method, associating 
the dispersive dielectric permittivity to a polarization current, 
as detailed in [12] and [26], is employed. For the dielectric- 
permittivity e(w), the Drude model is used, assuming an e/”? 
time-dependence 


|Q1Qal sin r = HRtt/? — pyntt/2, 3) 


Det! = pr + Hy.) (4) 
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—-— : Periodic boundaries 


Fig. 4. Implementation of the periodic boundary conditions in the triangular- 
mesh FDTD. 


B. Implementation of Periodic Boundary Conditions 


In a 2-D periodic structure of spatial period d, and dy 
along the x- and y-axes of a rectangular-coordinate system, 
phasor field components that are one period away in either 
direction differ only by a constant attenuation and phase-shift 
term exp(—jk,d,) and exp(—jk,d,), respectively, where k = 
ik, + Yky is a Bloch wave vector. This frequency-domain re- 
lationship is translated into the time-domain via the sine—cosine 
method of [20]. Note that the sine—cosine method enables the 
concurrent extraction of multiple resonant frequencies within 
the simulated frequency bandwidth that correspond to a Bloch 
wave vector enforced through the periodic boundary conditions 
[22]. As an excitation, a Gaussian point source, with a band- 
width covering the frequency range of interest is applied. For a 
fixed Bloch wave vector &, the time-domain waveforms of fields 
are sampled within the unit cell. Then, their Fourier transform 
(from the time-domain to the frequency-domain) reveals the 
resonances that represent the modal frequencies w(k). 

The following considerations are specifically pertinent to the 
implementation of periodic boundaries in the triangular FDTD 
mesh. Consider, for example, a wave propagating along the 
y-axis, with a Bloch wave vector k= gky, in a periodic struc- 
ture of the unit cell shown in Fig. 4. Note that although the 
magnetic-field node H, , is outside the unit cell, it is connected 
to the magnetic-field node H,,, at a distance d, along the 
y-axis inside the cell through the periodic boundary condition 


A,» = Aza exp(—jkydy). (6) 


Then, the tangential electric-field component along DF can 
be updated, according to the stencil of Fig. 3. However, the 
triangle enclosing the H,» node has to be the same as the 
one enclosing the H,, node for the two periodic bound- 
aries to be identical—not just physically but numerically as 
well. Therefore, the triangle DFG is generated by translating 
ABC by dy along the y-axis. Evidently, the segmentation of 
the periodic boundaries is also identical; if the left boundary 
(y = 0) is divided in segments P; P2,...,Py and the right 
boundary (y= d,) is divided in segments Q1Qz2,...,Qw, 
then |P; P11] = |Q:iQi+1|,2 = 1,2,...,N — 1. Similarly, the 
discretization of the lower and upper boundaries is identical. 
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Fig. 5. 


Equivalent circuit of the proposed triangular-mesh FDTD. 


Finally, modal-field distributions can be obtained by Fourier 
transforming sampled field values over a mesh of points. Note 
that the electric field at an arbitrary location within the compu- 
tational domain can be interpolated by using Whitney 1-forms 
[27], knowing the tangential component along each edge. 


C. Stability Analysis 


The triangular-mesh-based FDTD of this paper can be 
proven to be conditionally stable for a Drude medium by the 
equivalent-circuit approach of [28]. Therefore, (1) and (3) are 
considered at the sinusoidal steady state and discretized in 
space (with the triangular mesh employed in this paper) but not 
in time. The resulting system is called semi-discretized. The 
geometry and nodes involved are shown in Fig. 5. Magnetic- 
field components H,, and H,, are sampled at the centroids of 
triangles ABC and BCD, which are referred to as Q; and Qo, 
respectively. Also, the tangential electric-field components at 
the centers of CA, AB, and BC are denoted as E;,, F;,, and 
E,,. The spatially discretized sinusoidal steady-state version of 
(1) and (3), involving the aforementioned field components, is 


— jopAascH:, 

—E,, -CA+E,,:AB+E,, - BO (7) 
jwe(w)|Q1Qa| sin TEL, 

= H,, — Hz, (8) 


where Axygzc is the area of the triangle ABC. Introducing 
VY, = Et, CA, Vo = Et, AB, V3 = Ei, - BC, Va = Hz 
and V; = H,,, (7) and (8) can be written as 


1? 


—Ya4V4 = GyaVi + GoaV2 + Goa V3 (9) 


¥33V3 = G34V4 + G35 V5 (10) 


where Y44 = jwuAapc is the admittance of a capacitor, 
G4 = Go4 = G34 = 1 and G53 = —1 are gyrators, and Y33 = 
jwe(w)||Q1Q2|sin7/|BC]. If edge BC is within a lossless 
non-dispersive medium, Y33 is simply the admittance of a 
capacitor. If BC is within a lossless Drude medium, Y33 
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Fig. 6. Geometry of the SPC lattice of [9]. 


corresponds to a parallel LC’ circuit, while if it lies on the 
interface between two media, it can be characterized by an 
averaged €(w), which also corresponds to a parallel LC circuit. 
Therefore, the circuit analog of the semi-discretized system is a 
lossless passive network of gyrators, capacitors, and inductors. 
As aresult, the FDTD update equations stemming from the time 
integration of the semi-discretized system are not associated 
with late-time instability [28]. 

The formal stability condition for the triangular-mesh FDTD 
can be determined, via a standard methodology [18], to be 


2 


At < ———_. 
max; |w;| 


(11) 
In (11), w; are the eigenvalues of the semi-discretized system 
(the FDTD equations before the discretization of the partial 
derivatives with respect to time). The exact calculation of these 
eigenvalues is an unnecessary overhead for the computations 
of this paper. Instead of performing these, the time step is 
determined by the condition 


min; l; 


At=C 


(12) 
c 

where J; is the length of the ith edge of the triangular mesh, 
and C’ is an empirical constant: typically < 0.5. All simulation 
results that are presented next have been deduced with the con- 
stant C’ of (12) set to 0.1 to ensure that mesh nonuniformities, 
further reducing the distance between neighboring field nodes, 
would not render the scheme unstable. Indeed, stability has 
been achieved in all cases. Furthermore, no late-time instability 
has been observed for up to several millions of time steps, 
which is in agreement with the predictions of the analysis of 
this section. 


II. PERIODIC ANALYSIS OF THE SPC LATTICE 


A. Numerical Results: Rectangular- and 
Triangular-Mesh FDTD 


One of the motivating applications of this paper, namely the 
SPC lattice of [9], is considered next. The specifications of the 
geometry, as shown in Fig. 6, are normalized to the plasma 
frequency of the silver rods; the lattice period is a = c/wp, 
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Fig. 7. Power spectral density of H,, deduced by a rectangular-mesh-FDTD 
method for various sizes of Yee cells and Bloch wave vectors k = tke. 


Legends indicate ka values, where a is the lattice period. (a) a/100. 
(b) a/200. (c) a/400. 


and the radius of the cylinders is r = 0.45a. The SPC lattice 
was modeled in [9] by both quasi-static and full-wave (finite 
element) means. On the other hand, the triangular FDTD 
technique is employed here for the periodic analysis of the 
same structure. Recall that the FDTD process for this type of 
analysis proceeds by first fixing the phase shift between the 
periodic boundaries to correspond to a wave vector within the 
Brillouin zone and then sampling and Fourier transforming a 
field component inside the unit cell. The resonant frequencies 
associated with the wave vector at hand are determined by the 
position of field resonances in the frequency domain. 

To obtain more clear insights to the usefulness of the 
proposed formulation, results deduced via the conventional 
rectangular-mesh FDTD, for frequencies up to 0.7w,, are 
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Fig. 8. Power spectral density of Hz, deduced by the proposed FDTD 
technique for various minimum sizes of triangle edges and Bloch wave vectors 


k = &kz. Legends indicate kza values, where a is the lattice period. (a) a/20. 
(b) a/50. (c) a/100. 


presented first. In particular, Fig. 7 depicts the Fourier transform 
of H, sampled inside the unit cell of the SPC lattice for 
several Bloch wave vectors k = #k, (along ' — X). Although 
clear resonances can be seen up to 0.38 w,, higher frequency 
resonances present noise-like rather than resonant pattern. In 
addition, convergence of the resonant frequency peaks that 
appear above 0.38 w, cannot be achieved by mesh refinement. 
On the other hand, the triangular-mesh FDTD allows for a 
quickly convergent calculation of these resonant frequencies, 
even at a relatively coarse mesh, as shown in Fig. 8. 

The periodic analysis results derived by the triangular-mesh 
FDTD are summarized in Fig. 9, which compares the present 
analysis to the one of [9] (which provided the  — X part of 
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Fig.9. I — X part of the Brillouin diagram of the SPC lattice of [9]. Results 
from [9] are also shown for comparison. 


the Brillouin diagram of the SPC lattice). While the modes 
predicted by the analysis of [9] have also been found here, two 
additional ones have been extracted as well, revealing a richer 
band structure of the SPC lattice than previously thought. The 
nature of these modes is further illustrated by their magnetic- 
field magnitude patterns (all corresponding to k,a = 1, ky = 0) 
that are depicted in Fig. 10. The first mode (which is forward) 
shows a relatively uniform field distribution across the circular 
cross section of the plasmonic cylinder. The other four mode 
patterns are characterized by strong field localization at the 
air—cylinder interface, which is representative of their surface 
plasmon nature. 

Only the third of these modes is backward due to its slightly 
negative group velocity. This mode was used in [9] to realize a 
“perfect lens.” However, since this mode is closely surrounded 
by forward ones, the lensing effect cannot be maintained over 
a practically significant bandwidth. For completeness, the full 
Brillouin diagram of the SPC lattice is shown in Fig. 11. 


B. Discussion 


The comparison of the rectangular to the triangular FDTD 
results indicates that both can accurately determine the fun- 
damental forward mode of the SPC lattice. In fact, the clear 
Fourier transform peaks of Fig. 7, up to 0.38 w,, correspond 
to this mode. From that point on, the rectangular FDTD 
becomes poorly convergent due to the excitation of higher 
order plasmonic modes. The resonant frequencies of these are 
significantly perturbed by small errors related to staircasing 
approximations of the dielectric permittivity at the air—cylinder 
interface because of the strong surface confinement of their 
field patterns shown in Fig. 10, in contrast to the relatively 
uniform pattern of the fundamental mode. These perturbations 
do not converge with mesh refinement, as the latter simply 
redefines the air—dielectric interface modifying the associated 
surface plasmons and their frequencies found by FDTD. Such 
issues do not accompany the proposed approach, which can be 
a method of choice for studies involving surface plasmons at 
non-rectangular interfaces. 
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Fig. 10. Magnetic-field vector and magnitude plots for the SPC lattice modes 
of Fig. 9 obtained by the triangular-mesh FDTD. In all cases, kya = 1 and 
ky = 0. (a) First mode: w = 0.285 wp. (b) Second mode: w = 0.532 wy. 
(c) Third mode: w = 0.597 wy. (d) Fourth mode: w = 0.616 wp. (e) Fifth 
mode: w = 0.637 wp. 
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Fig. 11. Full Brillouin diagram of the SPC lattice of [9], as determined by the 
triangular-mesh-FDTD method. 
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Fig. 12. Isofrequency surfaces of the fundamental mode of the SPC lattice 
of [9]. The isosurface values correspond to the ratio w/w. 


IV. SPC LATTICE AS AN OPTICAL ANALOG OF 
MICROWAVE NRI-TL MEDIA 


Recently, a uniplanar NRI-TL medium inspired by the series 
of the Transmission-Line Matrix (TLM) technique [29] was 
presented in [23]. The authors of [23] also suggested an analogy 
between their microwave lattice and a triangular lattice of 
plasmonic nanospheres operating at optical frequencies. Since 
the 2-D version of this lattice is equivalent to the SPC lattice 
of [9], the periodic FDTD analysis of the latter can be used to 
evaluate this potential analogy. 

To this end, the isofrequency surfaces of the fundamental 
mode of the SPC lattice are plotted in Fig. 12. It is noted that the 
curvature of the isofrequency surfaces around the M-point sug- 
gests the propagation of a relatively broadband backward wave, 
which had not been pointed out before. While the formation of 
this backward wave is reminiscent of similar effects in photonic 
crystals [30], it is important to note that the unit cell of the 
periodic structure under study is deeply subwavelength (about 
A/16) at the frequency band where this wave is supported. It is 
also noted that the bandwidth of this wave is significantly larger 
than the SPC backward-wave mode that was discussed earlier. 

Finally, the electric-field pattern of this mode is shown in 
Fig. 13, which also indicates its circulating nature which is 
similar to the electric field supported by the cluster of plasmonic 
spheres studied in [6], as well as the uniplanar NRI-TL medium 
of [23]. This observation is also in agreement with the finite- 
element simulation results of [31]. 


V. CONCLUSION 


A simple, explicit, and stable formulation of the FDTD 
technique employing a triangular mesh has been shown to be 
a well-convergent efficient tool for the periodic analysis of 2-D 
arrays of plasmonic nanorods that are possibly useful for optical 
implementations of NRI media. Based on this technique, two 
contributions were made. First, the SPC lattice of [9] was 
fully characterized, and its rich plasmonic band structure that 
severely limits the bandwidth of a plasmonic backward-wave 
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Fig. 13. Electric-field vector and magnitude of the fundamental mode of the 
SPC lattice of [9]. Four cylinders are included to show the circulating pattern 
of the mode. 


band, which was previously employed for the realization of a 
“perfect lens,” was computed. Second, a relatively broadband 
backward wave, with a circulating electric-field pattern, was 
shown to be supported by the SPC lattice around the M-point 
of its Brillouin diagram. Therefore, this lattice may lend itself 
to optical NRI-TL implementations and associated negative- 
refraction applications. 
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